STUDY AREA

The study was done around Abisko, in northernmost Sweden (68°20′ N, 18°50′ E), on the north and the south side of the lake Torne Träsk (Fig. 1). The area is dominated by subalpine shrub heaths (e.g., Empetrum nigrum ssp. hermaphroditum, Betula nana, Arctostaphylos alpinus, A. uva-ursi, Vaccinium myrtillus, Phyllodoce caerulea, Carex bigelowii, Juncus trifidus) and open mountain birch forest (Betula pubescens ssp.). Mean altitude is 724 m a.s.l. (range 349–1365 m), with the treeline situated ≈650 m a.s.l. The southern side of the lake Torne Träsk receives 304 mm of precipitation yr⁻¹ (Alexandersson et al., 1991) and the north side of the lake twice as much (Sonesson and Hoogesteger, 1983). The mean annual air temperature in Abisko is −0.8° C, and the temperature of the warmest month, July, is 11.0° C (Alexandersson et al., 1991).

SATELLITE AND ANCILLARY DATA

We used an IRS 1D-LISS scene acquired 1 September 1998 (path 17, row 17). The image was orthocorrected to the Swedish National Grid (RT90) to within .5-pixel accuracy and resampled to 20 × 20-m pixels using cubic convolution (OM&M-Observation, Mapping and Monitoring AB, Stockholm). The IRS 1D is a high-resolution, multispectral image acquired from LISS-3 and mounted on an Indian satellite. The sensor collects radiation in 4 bands: Green (0.52–0.59 μm), Red (0.62–0.68 μm), Near Infrared (NIR) (0.77–0.86 μm), and Middle Infrared (MIR) (1.55–1.70 μm). The spatial resolution is 23 m for the first 3 bands and 70 m for the MIR band. Upon acquisition of the image, the sun elevation was 21 degrees above the horizon, and the sun azimuth was 163.8 degrees.

The elevation was derived from a digital elevation model (DEM) with 50-m resolution (National Land Survey) and resampled to 20 × 20 m using cubic convolution. Slope (0–90°) and aspect (0–360°) were extracted from a standard package (ERMapper 6.1). The curvature was extracted using a program in C that determines the second-order differentiation of the DEM. The first-order differentiation of the DEM gives the slope, and the second-order differentiation gives the curvature. The curvature can thus be interpreted as the variance in the slope, and this measure was found to be highly correlated with rock outcrops and microtopography (O. Hagner, pers. comm., August 2001). The curvature varies from 20 (concave) to 25 (convex) for the whole area.

FIELD DATA

Sample plots were laid out systematically on the northern and southern side of the lake Torne Träsk between 20 June and 2 September 1997. In total, 792 10-m-radius circular sample plots were positioned with a GPS (Trimble GeoExplorer ver. 2.11). The GPS data were differentially postprocessed to bring up the positional accuracy to <10 m. The sample plots were grouped into clusters of 9, forming a regular spaced grid of 60 m with clusters located 1 km apart (Fig. 1) (Dahlberg et al., 1998). In one north-south transect the clusters consisted of 15 plots, with distance between individual plots ranging from 30 to 500 m and distance between clusters of 500 m. Occurrence of reindeer on each sample plot was estimated by means of fecal pellet group counts, with 1 group consisting of ≥10 pellets. We used occurrence of readily identifiable pellet groups as a binary (0, 1) variable instead of using the actual number of pellet groups, reasoning that this approach provides a more stable, albeit conservative, measure of occurrence due to differences in detectability and decomposition rate across cover types. The field-sampled area on the northern side of the lake Torne Träsk is primarily summer and fall grazing ranges for reindeer, whereas the southern side is used primarily during spring and fall (County Administration, Norrbotten). The feces count thus represented the accumulated use of the area by the reindeer over the vegetation period.

The vegetation cover in each plot was classified according to the Swedish vegetation map (Liberkartor, 1981).

HABITAT SELECTION BY REINDEER

To derive the cover type most selected by the reindeer in terms of relative use (occurrence of fecal pellet groups), we employed Manly's selectivity index (Manly et al., 1993). This selectivity index can be described as the likelihood for a particular cover type to be selected, given that all cover types are equally available. Manly's selectivity index has the following form:

STATISTICAL MODELING

Elaborating the Final Model

Once the best of the 3 models was estimated, we examined the autocorrelation in the response variable. A variogram of the response variable was drawn prior to its being converted into a binary variable (Fig. 2). As can be seen from the variograms, the correlation between the response variable and the location of the field data plots is fairly high for distances up to 2 km. Given that our plots are much closer together than this, we decided to include a measure of the spatial correlation in the logistic model.

When a classical regression model is applied to data that exhibit autocorrelation, the estimation of the parameters is not efficient, and the estimates of the variances are biased. If the autocorrelation is positive, the mean square error is underestimated, resulting in an optimistic measure of fit: R² is higher than the true value. Furthermore, the standard errors are underestimated, which will lead to tests of significance that are misleading. The reverse occurs when the data have a negative autocorrelation (Upton and Fingleton, 1985; Long, 1998; Lark, 2000). The same problem plagues logistic regression. This problem has been addressed by a number of authors who have proposed the autologistic model (Besag, 1974; Ripley, 1981; Cressie, 1993; Augustin et al., 1996). In this model, the spatial autocorrelation is represented by an autocovariate, which is a weighted sum of the neighboring pixel values. The weights selected here are the reciprocal of the Euclidean distance between two cells.

where

where x_i represent the explanatory variables, f_i are appropriate transformations, α and β are the parameters that need to be estimated, k_i is a neighborhood around cell i, and d_ij is the Euclidean distance between cell i and cell j.

As suggested by Augustin et al. (1996), one can estimate the autocovariate through an iterative process. A simple logistic regression is modeled using the data, and the probability of occurrence of particular cover types is estimated for all cells where no field data exist. The autocovariate is then estimated, and an autologistic model is fitted with the autocovariate as an additional explanatory variable. This process is iterated until it converges; in our case this meant convergence in the area under the ROC curve.

HABITAT SELECTION BY REINDEER

Occurrence of reindeer pellet groups (Table 1) showed that reindeer used snowbed vegetation most. Therefore, we selected snowbed vegetation for the modeling work. Snowbed vegetation occurred in 127 of the 792 plots (16%), but because snowbed vegetation occurred only in plots >700 m a.s.l., we excluded all plots below this altitude, leaving 398 plots for the analysis. Thus, our training data set comprised 300 plots and the evaluation data set 98 plots.

EXPLORATORY MODEL ASSESSMENT

Because a binormal model (Pearce and Ferrier, 2000) did not provide a good fit to our data, we derived the ROC curve by fitting a polynomial of order 5 through the data. The area under the curve (AUC) was then obtained by integrating the polynomial function. The AUC for model 1 and model 2 was 0.727 and 0.688, respectively. For model 3, i.e., the combined model, the AUC was 0.756. Model 3 combined with the autocovariate (model 4) produced an AUC of 0.771. The ROC curves for all 4 models are overlaid in Figure 3. One observes that the true positive fraction increases steeply for low cut-off probabilities, slightly more so for model 4. The corresponding cross-validated accuracy rates (at a cut-off probability of 0.25) were 67.3, 55.0, 74.5, and 74.5%, respectively. This cut-off probability was determined based on the prior distribution of snowbed vegetation in the test area (Hosmer and Lemeshow, 1989).

PARAMETER ESTIMATES AND MODEL EVALUATION

The estimates of the parameters of model 3 and model 4 are listed in Table 2. Note that an increase in Band 1, Band 3, and Band 4 will increase the expected probability of snowbed vegetation. Similarly, the coefficient of NDVI is negative. The parameter estimate for curvature is also negative, which means snowbed vegetation is more likely in a concave landscape. Note that no interaction terms were significant. The contribution of each variable to the overall reduction of the deviance for each model is given in Table 3. It is clear that the most important variables are Band 1, Band 3, NDVI, and curvature. The signs of the variable have not changed, but some slight changes occur in most of the values of the parameters; some have increased and others have decreased. These changes are related to how well the different variables are affected by spatial information. The autocovariate term contributes least to the deviance decrease, but it is significant at the 95% level. The probability surface map corresponding to the autologistic model is given in Figure 4.

MODEL EVALUATION

There exist several approaches for modeling data that are spatially autocorrelated. On the one hand we have regression types techniques where the spatial autocorrelation can be implicitly (Upton and Fingleton, 1985) or explicitly (Mardia and Marshall, 1984) modeled. On the other hand we have geostatistical models, which originated from the mining industry (Cressie, 1993). However, before one can choose among such a panoply of statistical tools, one has to make a decision about the level at which the data can be modeled. One could consider our data to have 2 hierarchical levels. At the upper level, we have the original data, with the proportion of snowbed vegetation in the plot. This proportion varied in parts of 10 from 0 to 100. At the next level, we have an indicator variable that gives the presence and absence of snowbed vegetation. Given the amount of data available and the complexity of the problem, and after preliminary analysis, we found that we could model only the indicator variable with an acceptable level of accuracy and stability. Given this decision, the geostatistical approaches became less attractive since the theory of indicator kriging is not very well developed (Journel, 1983; Jiménez-Espinosa and Chica-Olmo, 1999). Furthermore, the autologistic approach provides an output that is a probability, and the algorithms for estimating such models are fairly well developed and tested.

The logistic statistical model is a very useful approach for incorporating data of different types from several sources. The logistic model was found to be similar to the linear discriminant function (Maddala, 1993). Manel et al. (1999) found that logistic regression outperformed discriminant analysis and artificial neutral networks when applied to independent data sets. Given the close proximity of the plots within a cluster and the high correlation between the neighbors, the autologistic model is useful and performed fairly well. Including terrain characteristics and ratios of the bands greatly improved the efficiency of the models. Furthermore, because the logistic model provides probability of occurrence of snowbed, it acts as a soft classifier, the results of which can be further improved by some appropriate postprocessing algorithms (Foody, 2000). Given the algebraic consequence that space is represented as an additive covariate, the improvement of the autologistic over the simple logistic model is mostly observed as a smoothing of the probability image. This can be observed when transects along the two probability maps are compared. In Figure 5, we see a difference image for a subset of the image around the test site, produced by subtracting the probability maps produced by model 3 and model 4. It should be noted that the probabilities estimated from model 4 (autologistic model) are always lower than those produced by model 3 because of the averaging effect of the autocovariate.

CONCLUSIONS

We developed a probability surface algorithm that proved useful for mapping a cover type having small areal extent and/or a scattered distribution. In the future, we hope to examine how pixel-wise probability classification maps of vegetation cover types can be integrated with maps based on coarse-scale classifications of mountain vegetation so as to maximize the information content. This integration can be particularly useful because coarse-scale classifications exist for large areas and contain a wealth of information covering the ecoregion (dominant cover types) (e.g., He et al., 1998). Pixel-wise classification allows fine-scale assessment of distribution of important vegetation cover types. Combined with coarse-scale classifications, pixel-wise classification schemes could be useful for assessing the distribution and abundance of limiting habitat resources, e.g., for reindeer management, and environmental monitoring. Because probability maps do not depict the exact occurrence of particular cover types but rather their relative occurrence, they are best suited for landscape-level assessments (10³–10⁴ ha). Deriving fine-scaled probability map surfaces requires some amount of field sampling, which may be difficult to achieve, particularly in remote parts of the Arctic. On the other hand, to be taken seriously, remote-sensing assessments must be based on solid ground truth data.